Two-band second moment model and an interatomic potential for caesium 
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A semi-empirical formalism is presented for deriving interatomic potentials for materials such as 
caesium or cerium which exhibit volume collapse phase transitions. It is based on the Finnis-Sinclair 
second moment tight binding approach, but incorporates two independent bands on each atom. The 
potential is cast in a form suitable for large-scale molecular dynamics, the computational cost being 
the evaluation of short ranged pair potentials. Parameters for a model potential for caesium are 
derived and tested. 

PACS numbers: 34.20.Cf , 61.72.-y, 64.60.-i 
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I. INTRODUCTION 

Semiempirical models for metallic binding have had a 
long and successful history in computer modelling. The 
most significant development came in the mid eighties 
with the implementation of 'embedded atom' potentials 
jjj based loosely on density functional theory [|| and 'N- 
body' potentials || based on the tight binding second 
moment approximation pi. 

Although the rationale for these potentials suggests ap- 
plicability to free-electron and transition metal systems 
respectively, the implied functional form is similar in each 
case. A large number of successful parameterisations of 
the functional forms suggested by these methods have 
been used, and they have established themselves as the 
standard method for modelling metallic systems. 

In the second moment approximation to tight binding, 
the cohesive energy is proportional to the square root 
of the bandwidth, which can be approximated as a sum 
of pairwise potentials representing squared hopping inte- 
grals. Assuming atomic charge neutrality, this argument 
can be extended to all band occupancies and shapes 
For simplicity, consider a rectangular c£-band of full width 
W centred on the atomic energy level Eq . The energy for 
this band relative to the free atom (the bond energy) is 
given by: 
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E—dE= n(n-N) (1) 



where n is the occupation of the band and N the capacity 
(for d-bands N=10, for s-bands N=2). The band width 
is defined as the square root of the second moment of the 
density of states. If we project the density of states onto 
each atom, and assume that each atom is charge neutral 
then the energy becomes: 
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and the second moment of the density of states can be 
calculated as the sum of the squares of the hopping in- 
tegrals to nearest neighbours. This latter operation can 
be written as a sum of pair potentials, 
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(3) 



The enormous success of Finnis-Sinclair (and embedded 
atom-type) potentials arises from their extreme compu- 
tational efficiency: essentially they are no more compu- 
tationally expensive than a conventional pair potential. 
This allows them to be applied to extremely large scale 
simulations, addressing complex geometries and phenom- 
ena which are intractability with more accurate quantum 
mechanical models. 

The computational simplicity follows from the formal 
division of the energy into a sum of energies per atom, 
which can in turn be evaluated locally. Here we follow a 
similar philosophy, seeking to incorporate as much of the 
relevant physics as possible, without increasing compu- 
tational cost. 

The alkali and alkaline earth metals appear at first 
glance to be close packed metals, forming fee, hep or 
bec structures at ambient pressures. However, compared 
with transition metals they are easily compressible, and 
at high pressures adopt more complex "open" structures 
(with smaller interatomic distances). The simple picture 
of the physics here is of a transfer of electrons from an s 
to a d band || |9| , the d-band being more compact but 
higher in energy. Hence, at the price of increasing their 
energy, U, atoms can reduce their volume, V, and since 
the stable structure at OK is determined by minimum en- 
thalpy, H=U+PV at high pressure this transfer becomes 
energetically favourable. The net result is a metal-metal 
phase transformation characterised by a large reduction 
in volume (and often also in conductivity), the crystal 
structure itself is not the primary order parameter, and 
in some cases the transition may even be isostructural. 

The mechanism of the phase transition is unknown. 
Although it is obvious that an isostructural phase tran- 
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sition is accompanied by an instability of the bulk mod- 
ulus, the onset of this instability may occur after that 
in the shear modulus. Thus the mechanism may involve 
shearing rather than isostructural collapse, particularly 
if a continuous interface between the two phases exists, 
as in a Shockwave 11(1 . 



II. THE TWO BAND MODEL 

A. Energies 

The Finnis Sinclair formalism has been successfully 
used to study a large variety of systems. However, we 
are interested in systems in which electrons change from 
one type of orbit to another and in particular, from an 
s-type orbital to a d-type orbital as the sample is pres- 
surised. 

Consider therefore two rectangular bands of widths W\ 
and W 2 as shown in figure [l] with widths evaluated using 
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N/W 2 + N/W z 
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equation |. The bond energy of an atom may be written 
as the sum of the bond energies of the two bands on that 
atom as in equation |^, and a third term giving the energy 
of promotion from band 1 to band 2 (see equation^): 

EWa W i2 
T^nii(nii-N 1 ) + ——n l 2(n l 2-N 2 )+E prom 
ZlXi 2iv2 

i 

(4) 

where Ni and N 2 are the capacities of the bands and nn 
and rii2 are the numbers of electrons in each band on the 
i th atom. 

For an element, enforcing charge neutrality, the sum of 
the numbers of electrons is necessarily constrained to be 
equal to the total number of electrons on an atom, 



n-n + n l2 = T . 



(5) 



The difference between the energies of the band centres 
«i and a 2 is assumed to be fixed. The values of a cor- 
respond to the appropriate energy levels in the isolated 
atom. Thus, a 2 — u\ is the excitation energy from one 
level to another. For alkali and alkaline earth metals, 
the free atom occupies only s-orbitals, so the promotion 
energy term is therefore simply 



n 2 (a 2 - a\) = n 2 E , 



(G) 



a 2 -w 2 /2 a 2 a 2 +w 2 /2 



where Eq — a 2 — a±. 

Thus the band energy can be written as a function 
of nil, na and the bandwidths (evaluated at each atom 
as a sum of pair potentials, within the second moment 
approximation) . Defining: 



FIG. 1: Schematic picture of DOS in rectangular two band 
model. Shaded region shows those energy states actually oc- 
cupied. 



Vi = na - n i2 
and using equation ^ we can write: 



(7) 



I 



Although this expression looks unwieldy, it is in fact com- 
putationally efficient, requiring for its evaluation only two 
sums of pair potentials for W (see equation |^) and a 
minimisation at each site independently with respect to 
r/i. Since they are local variables, it is possible to write 
closed form expressions for r\i which minimise the total 
energy (see equation |l4|) . 

In addition to the bonding term, the Finnis-Sinclair 
form includes a pairwise repulsion between the ions, 
which is due primarily to the screened ionic charge and 
orthogonalisation of the valence electrons. In the present 



case, this pair potential should be a function of r/i. In 
keeping with maintaining locality of the energy we write 
the pairwise contribution to the energy in the intuitive 
form, as the sum of two terms, one from each "band" , 
proportional to the number of electrons in that band 



V{r i3 ) = (n a + nji)Vi{r i: j ) + (n l2 + n j2 )V 2 {r i:j ) (9) 



We rearrange this to give the energy as a sum over atoms 
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Tin ^ Vi(rij) + n l2 ^ v ^( r ij) 



This depends on rji, but as described above the rji take 
the values which minimise the energy and one can solve: 



(10) 



dUt, 



The total energy is now simply 

Utot — Up a i r H~ Ubond • 



= 



(12) 



(11) explicitly for rji independently at each atom, whence 



Via 



NiN 2 



W a N 2 + W a N! 



II /J. - I' i-l — T ( -J^- ] + 'lEu — "21 ,!./„//, + -l : ,2.p<,ir 



(13) 



Depending on the number of electrons in the system 
it may not be possible to realise this. Charge neutrality 
requires that \rji\ cannot be greater than the total number 
of electrons T per atom. The fixed capacities of the bands 
(JVi and N 2 ) can also prohibit the realisation of T)i . It 
is therefore necessary to limit the values which rji may 
have: 



V 



min(T, 2Ni 
max(— T, T - 
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T), if mo > min(T, 2Ni ~ T) 
2N 2 ), if r)i < max(-T, T - 2N 2 ) 
otherwise 

_ (14) 
where r) io is given by equation p~3l The expressions for 

rji involves only constants and sums of pair potentials, 
and can be evaluated independently at each atom at a 
similar computational cost to a standard many-body type 
potential. 



B. Forces 

Since the energy is variational in rji within this model, 
we have 



dU tt 



drjt 



= 



(15) 



and the force on the i atom is 
dUtot 



dvi 
dU tot 



dri 

dUtot 



drj dri 



(16) 



Hence the force is simply the derivative of the energy at 
fixed rj. Basically, this is the Hellmann-Feynman theorem 
p| which arises here because rj is essentially a single 
parameter representation of the electronic structure. 



This result means that, like the energy, the force can 
be evaluated by summing pairwise potentials. Hence the 
model is well suited for large scale molecular dynamics. 

The force derivation is somewhat tedious (see ap- 
pendix |b|), the result being: 

f i = + WilMN) + {Wi2 + W ]2 )4>' 2 {r l3 ) 



(riii + nji)V'i(rij) - (n l2 + n j2 )V' 2 (nj) 



(17) 



with 



±N b Wib 



and where the numbers of electrons in each band, na, 
have been calculated analytically using equations [l3|, [l4] 
and 0. The contribution to the force from each neighbour 
acts along the direction of the vector between the atoms. 



III. A MODEL POTENTIAL FOR CAESIUM 

One application of the model is to simulate the 
pressure-induced isostructural phase transition in mate- 
rials such as caesium. Here the transformation arises 
from electronic transition from the s to the c?-band. Cae- 
sium adopts the bec structure (Cs I) at ambient pres- 
sure, transforming under slight additional pressure to fee 
(Cs II) and under further increase transforming isostruc- 
turally to fee (Cs III) Here we seek to represent 

the volume collapse Cs II — ► Cs III transition. 



A. Parametrisation 

To make a usable potential, the functional forms of 
(j> and V must be chosen. Although this is somewhat 
arbitrary, the physical picture of hopping integral and 
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screened ion-ion potential suggests that both should be 
short ranged, continuous and reasonably smooth. 

In the present case, we are interested in the gross fea- 
tures of the model, so we adopt simple forms fitted to the 
cohesive energies and volume of the high and low pressure 
polytypes of caesium. For more complex applications, it 
may be necessary to fit other properties. 

Following Finnis and Sinclair, we choose for the hop- 
ping integral 



M r ij) 




if r»j < d b 
if m > db 



and for the pairwise part, following Lennard- Jones 



(18) 



(19) 



where b labels the band. 

We select the promotion energy Eq to be that required 
to promote an electron from the 6s level into the bd level 
of an isolated atom |l4|, The band capacities are 

N s = 2, Nd = 10 and the total number of electrons per 
atom is T — 1. The cutoff radii d s and dd are chosen to 
be between the second and third nearest neighbours and 
between the first and second nearest neighbours respec- 
tively. The remaining four parameters (A s , Ad, C s , Cd) 
are available for fitting. 

Ab initio calculations using pseudopotentials, plane 
waves and the generalized gradient approximation |l6| , 
|l7| show that the energy-volume relations for bcc and 
fee caesium are almost degenerate, so we have fitted 
the remaining parameters to the atomic volumes of Cs I 
(115.9 A 3 /atom), Cs II (67.5 A 3 /atom) and Cs III (48.7 
A 3 /atom) |l4j as well as the transition pressure between 
phases II and III (4.3 GPa) and the cohesive energy 
of phase I (-0.704 eV). The cohesive energy is the sum of 
the heat of formation and the heat of vaporisation [Q . 

Initial estimates for the parameters were determined 
using a symbolic mathematics package. The final pa- 
rameters were arrived at by an iterative process: Least 
squares, conjugate gradients and ab oculo minimisation 
techniques were used to determine the best parameters 
for particular cutoff radii. The cutoff radii were then 
adjusted by hand to improve the fit. This process was 
repeated until the optimum fit was achieved. The final 
parameters are given in table Q. 

Figure [2] shows the energy-volume curves for the fee 
and bcc structures calculated using the model. The fee 
structure is stable everywhere compared to the bcc struc- 
ture. Experimentally however, Cs I has bcc structure 
with an equilibrium volume of 116 A 3 . The present po- 
tential does not have a low pressure bcc phase. However, 
at ambient pressure the bcc and fee curves are almost 
degenerate (0.02 eV/atom difference) in agreement with 
the ab initio calculations. The predicted equilibrium vol- 
ume of the ambient pressure phase agrees well with the 
experimental equilibrium volume of phase I of caesium. 



s band 


d band 


C a 


0.05617 eV 2 A~ 3 


c d 


0.1681 eV 2 A~ 3 


d s 


9.5097 A 


d d 


6.9189 A 


A 3 


2.4017 x 10 7 eV A 12 


A d 


3.7668xl0 6 eV A 12 


E 


1.19 eV 



TABLE I: Parameters for the two band model, obtained by 
fitting to phases II and III of caesium. 



Although the predicted volumes of phases II and III are 
larger than the experimental values, the predicted transi- 
tion pressure and volume collapse are in good agreement 
with experiment. 
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FIG. 2: Top: Variation of r/i with compression, showing s — » d 
transfer in model caesium (T — 1). Bottom: An energy- 
volume curve for the two band potential. The minimum lies at 
-1.3163 eV and 115.2 A 3 /atom. The gradient of the straight 
dash-dotted line is the experimental fcc-fcc transition pres- 
sure. In reality caesium also has a bcc-fcc phase transition at 
2.3 GPa, however first principles calculations show that these 
two structures are almost degenerate in energy at OK. The 
isostructural transition, which involves electron localisation 
and hence is not properly reproduced by standard electronic 
structure calculation, occurs at 4.3 GPa. 



B. Elasticity 

The elastic moduli are not fitted explicitly, so their 
behaviour represents a sensitive test of the model. Ap- 
plication of the pressure generalisation of the Born sta- 
bilitycriteria is rather confused in the present literature 
|i~9| ptif so we lay this out in appendix [a|. 

Since caesium exhibits an isostructural phase transi- 
tion the bulk modulus must formally become negative at 
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some volume (where the structure is unstable). It is less 
clear whether the other Born stability criteria will be 
violated: ab initio density functional perturbation the- 
ory calculations Jl7f| find an instability in the long range 
acoustic phonons: at zero pressure this is equivalent to 
an instability in C = (Cn — C\2)/2. 

With the present potential the analytic expressions for 
the elastic constants are complicated: there is no simpli- 
fication akin to equation [l5| for the second derivatives. 
Consequently, we evaluate the elastic constants numeri- 
cally from finite strain calculations, see figure | We find 
that the volume collapse is announced by a slight soft- 
ening of the bulk modulus which then goes negative in 
the unstable region. Although the shear and tetragonal 
shear decrease in the unstable region, neither actually 
goes negative. 
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40 50 60 70 80 90 100 110 120 
volume per atom (A 3 ) 

FIG. 3: Graph of elastic stiffness constants against volume. 
Mechanical instabilities occur for negative values of B44, 
- B 12 )/2 and (B u + 2Bi 2 )/3. Total energy is shown 
on the same figure but with respect to the right hand axis. 
The right most vertical line indicates the volume at ambient 
pressure. The other two vertical lines delineate the volumes 
in which the structure is unstable with respect to decompo- 
sition into coexisting phases whilst the straight dash-dotted 
line is the transition pressure from phase II to phase III. 



C. Defects 

One of the major features of many-body and embedded 
atom type potentials is their ability to describe defects 
such as surfaces, vacancies and interstitials without the 
constraints of pair potentials, for which the undercoor- 
dination defect energy is simply the sum of broken bond 
energies, less a small amount from relaxation of atomic 
positions. Thus the vacancy formation energy is typically 
the same as the cohesive energy, whereas in real materi- 
als it is typically less than half. The present potential is 
not fitted to any defect configuration, so it is of interest 
to see what it predicts - moreover the variety of envi- 
ronments associated with defects provides a good check 
against pathologies. 



surface 


energy me V / A^ 


relaxation lA 


relaxation 2A 


111 


5.938 


-0.042 


0.000 


001 


6.54 


-0.046 


-0.007 


011 


7.195 


0.024 


-0.015 


211 


6.952 


-0.122 


-0.045 



TABLE II: Surface energies, and atomic relaxations of the top 
two layers in angstroms 



configuration 


energy eV 


111 dumbbell 


2.178 


001 dumbbell 


1.776 


011 dumbbell 


1.975 


011 crowdion 


1.975 


vacancy 


0.545 



TABLE III: Interstitial and vacancy formation energies. The 
vacancy formation volume is 0.846Vo- 



The relaxed surface energies for four low-energy sur- 
faces are given in table [n]. As is usual with many-body 
potentials the energy is much lower than would be ex- 
pected from simple bond counting. For the (111) surface 
the energy for each atom in the plane is increased by 
about 1/9 of its cohesive energy, while each atom has 
1/4 of its bonds broken. There is also transfer of all 
electrons from ri-like to s-like states at the free surface. 
This also leads to an unusual outwards relaxation of the 
surface atoms. 

Like wise the vacancy formation energy and volume (ta- 
bic III) are rather typical of many-body potential results. 



The interstitial formation energy is especially low (ta- 
ble HI), around 1.8 eV depending on the orientation and 



the amount of relaxation allowed locally, meaning that 
thermal interstitial formation is possible. This low value 
is due to the transfer of electrons from s- to ri-band on the 
interstitial atom or dumbbell, similar to the high pres- 
sure behaviour. Consistent with recent ab initio calcula- 
tions in various elements [^l[ |2^ , ^3| , the calculated in- 
teratomic spacing of the dumbbell atoms is much smaller 
than in the bulk, (around 15%) and much smaller than 
is typical of standard EAM-type potentials. As a conse- 
quence of this, the associated strain fields are consider- 
ably smaller. 

For a detailed study of point defects in caesium, it 
would be appropriate to reparametrise the potential with 
point defects included in the fitting, but the good re- 
sults obtained here without such fitting suggest that the 
present model contains the right physics. 
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IV. EXTRAPOLATION TO TRANSITION 
METALS 

In principle, the current formalism should be appli- 
cable to ci-band metals. We do not intend to refit the 
potentials here, but by applying the parameters fitted 
for caesium with appropriate scaling for ionic charge, 
and simply varying the total number of electrons we re- 
cover the parabolic behaviour of the cohesive energy and 
bulk modulus which characterises the transition metal 
series. Considering that there are no fitting parameters, 
the agreement with experimental results is extraordinar- 
ily good. Moreover, the volume collapse phase transition 
exists only for N=l and N=2 (consistent with experi- 
ment). However, since no information about band shape 
is included, the sequence of crystal structures cannot be 
reproduced, so we consider here only the fee structures, 
calculating their "experimental" values from the experi- 
mental density. 

While the extrapolated potentials do not represent the 
optimal parametrisation for specific transition metals, 
the recovery of the trends across the group lends weight 
to the idea that the two-band model correctly reproduces 
the physics of this series. 




t i i i i i i i i i i i_ 

1 23456789 10 11 
number of sd electrons 



FIG. 4: Extrapolation of model to various 6s5d-band materi- 
als. The parameters Eo, A 3 , Ad, C 3 , Cd, d 3 and dd are taken 
from the fit to caesium, with the lengths being scaled accord- 
ing to the Fermi vector (Z _1//3 ) and the energies by (Z 1 ^ 2 ). 
The lattice parameters are shown by the solid line (calculated) 
and circles (experiment) and the cohesive energies by dashed 
line (calculated) and squares (experiment). Calculations are 
actually done at integer values of T and lines join these points. 



V. CONCLUSIONS 

We have presented a model to describe s — > d transfer 
within the framework of the empirical second-moment 



tight binding model. With a very simple parametrisa- 
tion, our model describes the isostructural phase tran- 
sition and associated s — > d transfer in Cs and allows 
study of the elastic instabilities which occur under pres- 
sure and uniaxial stress. Throughout, our formalism is 
guided by the principle of retaining as much electronic 
detail as possible within the constraint of the computa- 
tional complexity of short-ranged pair potentials 

Our potential form involves evaluation of sums of pair 
potentials, and minimisation of the energy with respect 
to band occupation. In the approximation of local den- 
sity of states evaluated at each atom, it becomes possible 
to formally write the energy of each atom in terms of the 
local band occupation, a single parameter rji. This vari- 
ational parameter can be explicitly eliminated from the 
energy and force expressions, meaning that energies and 
forces can be evaluated by summing pair potentials. 

To our knowledge, the present potential represents the 
most sophisticated energy function which has first deriva- 
tives which are analytic pair potentials, and the first ex- 
ploitation of the Hellman-Feynman theorem in empirical 
potentials to achieve this. It is thus uniquely well suited 
for exploitation in standard molecular dynamics codes. 

Extension of the formalism to transition metals is 
straightforward, and the caesium parametrisation gives 
a surprisingly good description across the group. While 
providing impressive proof of concept, this may not be 
of enormous practical use, however, since once there is 
a significant partial occupation of the d-band under all 
conditions the model is similar to the standard many- 
body potentials which are known to do a good job in 
describing transition metals. 

The precise stable crystal symmetry depends on the 
shape of the band structure, and is not therefore cor- 
rectly described in this rectangular band model. Many- 
body potentials never contain the correct physics to de- 
scribe phase transitions, e xcep t in the martensitic case 
of freezing in soft phonons p4 fmI - However, the crystal 
structure which has the lowest energy can be determined 
by judicious choice of functions V and (f>. Here the simple 
use of power law repulsion and cubic attraction leads to 
close packed structures. 

The two-band formalism can be easily combined with 
alternate parameterisations to get effective band shapes, 
and could be applied to ferromagnetic systems where the 
two bands represent separate spins, and the atomic Eq 
term is replaced by a term which favours maximum spin. 
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I 

APPENDIX A: ELASTIC INSTABILITIES UNDER PRESSURE - ENERGY DEFINITION 



There is a degree of confusion in the literature regarding the definition of elastic constants under pressure. To some 
extent this arises because of various definitions of strain (Lagrangian, Eulerian, volume conserving). Here we lay out 
the definitions in terms of energy. One important point to note is that under pressure the moduli which correspond 
to long wavelength phonons are volume conserving, while those corresponding to crystal stability are not. At finite 
pressure, this means they are different. 

At pressure, the important quantity for stability in energetic terms is the Gibbs free energy G. 

G = U -TS + PV (Al) 

where P is an externally applied hydrostatic pressure. U — TS is the Hclmholtz free energy F. The change in the 
Gibbs energy due to a distortion of the crystal is 

AG = AF + AP.V + PAV . (A2) 

If the hydrostatic pressure is applied by an external mechanism then AP is zero. Therefore 

AG = AF + PAV . (A3) 

Any arbitrary change in the unit cell be expressed in terms of the strain tensor e . Using Voigt notation, the matrix 
representation of the strain tensor is written as 

/ ex e 6 /2 e 5 /2\ 
e= e 6 /2 e 2 e 4 /2 . (A4) 
\e 5 /2 e 4 /2 e 3 / 

If we write the equilibrium cell as a matrix comprising the three lattice vectors Va= *2 r 3)> then any arbitrary 
strain give a new unit cell: 

V=Vo (1+ £ ) ■ (A5) 
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If the volume of the original cell is Vb = | Vb | , then the volume of the new cell is 

V/V = | y |/| Vo | (A6) 
= 1 + e\ + e 2 + e 3 + e x e 2 + e 2 e 3 + e 3 ei 
-e\l '4 - el/ '4 - 

+eie 2 e 3 - eief/4 - e 2 e\/A - e 3 eg/4 + e 4 e 5 e 6 /4 . 
Which can be expressed concisely in standard strain notation (i, j, k, I representing cartesian directions) as: 

-X- = e-u + - (2SijSki - S ik Sji - SuSjk) e tj e kl + 0(e 3 ) (A7) 

where Sij is the Kroneker delta and the usual implicit sum convention for tensors applies. The change in the Gibbs 
free energy may then be written 

PV 

AG = AF + Pea + — (26 l3 6 kl - 5 tk 5ji - 5 a 5 jk ) e tJ e M . (A8) 

We define Cy as the second order derivatives of the Helmholtz free energy with respect to the Voigt strains, and 
Bij as the second order derivatives of the Gibbs free energy with respect to the Voigt strains. The limiting case of long 
wavelength phonons is related to Cy which can be calculated from the dynamical matrix p^ ], the crystal stability 
criteria are related to the In the standard notation: 

P 

Bij k i — Cij k i + — (2Sij5 k i — S ik 8ji — SuSjk) ■ (A9) 

In Voigt notation: 

!C U - P/2 if i = j and i > 3 

Cij +P if % ^ j and i,j < 4 (A10) 

Cij otherwise 

The Born elastic stability criteria restrict what values various moduli may have if the crystal is to be stable p7| . In 
the case of external pressure the relevant free energy is G and the relevant moduli are B^ . Due to the Voigt symmetry, 
the elastic constant tensor has in general, only 21 independent components. In cubic crystals however this number is 
further reduced by symmetry to 3, namely Bn, B\ 2 and -B44. Consequently, there are three common stability criteria 
which impose lower bounds on the bulk, shear and tetragonal shear moduli of cubic crystals. These may be written 
in Voigt notation: 

i(Bn+2B 12 ) >0, 5 44 >0 and i(B u - B 12 ) > . (All) 
Whence the generalised Born stability criteria for the elastic constants of cubic crystals at pressure become: 

~(C n +2C 12 + 2P) > 0, 
C 44 -| > 0, 

\{Cxx-C 12 -P) > 0, (A12) 



where 



B n = C n , B 12 = C 12 + P and B 4A = C u - | . (A13) 



Calculation of bulk, shear and tetragonal shear moduli from energies 



In this paper we evaluate the elastic moduli by applying finite strain and measuring the change in total energy 
Utot- The finite strain method has the advantage over analytic differentiation of automatically compensating for 
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non-isotropic movement of the atoms. At zero temperature Utot is the same as the classical Helmholtz free energy F. 
At finite pressure, however, the system is not at a minimum of U to t with respect to strain. 

When applying strain to a system under pressure, the first order term in the Gibbs free energy vanishes: 

dG 1 d 2 G 1 d 2 G 

AG = ey h - en e kl - — - — = - ey e kl - — - — (A14) 

deij 2 ddijdeu 2 de^deu 

But the first order term for Helmholtz does not. For a cubic crystal, one can find (C\\ + 2Ci2)/3, (Cn — Ci2)/2 and 
C44 directly at any volume by simple distortions and measurement of F. 

For bulk modulus we apply a pair of Voigt strains, with e\ = e 2 = = e, the other elements being zero. Assuming 
we make a pair of small distortions about a volume V, the difference 

G(e) + G(-e) - 2G(0) = F (Bn + 2gia) 9e 2 + 0(e 3 ) . (A15) 

equivalently, using the Helmholtz free energy 

(C 11 + 2C ia ) = {F(e) + F(-e)-2F m 
3 9Ve 2 

For (Cn — Ci2)/2, the relevant distortion, volume conserving to first order, is e = e± = —e 2 and 

(di - C 12 ) (F(e) + F(-e) - 2F(0)) 



2 We 2 
For (G44, the relevant distortion is e = ee and 



0{e 6 ) . (A17) 



W.) + JT( e )-2F(0)) 



APPENDIX B: DERIVATION OF THE FORCE ON AN ATOM 



Since the cohesive energy is variational with respect to the electron distribution 77, the derivative of the energy 
with respect to the electron distribution is zero and so the force may be written as a partial derivative akin to the 



Hellmann-Feynman theorem JL1[ . From equation 16 



dU t , 



dr,. 



E 



9Ui 
dvi 



(Bl) 



where the sum over j includes j — i. From equations ^, ^ and ^, the cohesive energy of the j th atom in terms of the 
number of electrons in each band (riji and rij 2 ) is 



^rvm' 1 Wjl + ^\N- 2 



1 j W j2 + riji ^2 Vi(r jk ) + n j2 ^ V 2 (r jk ) + n j2 E 



and the force on atom i is thus 

U = - 



\- \njifnji _ A dWji n j2 ( n j2 _ \ dW j2 
4^ [ 2 V Ni J dvi 2 \N 2 ) dr, 



d \ ( d 

riji + nu—Vi(rij) J + \ nj 2 + n i2 — V 2 (r i:j ■) 



(B2) 



(B3) 



The width of band b on atom i is Wjb — \/ X^fe^j ^bi^jk) & n d its derivative is 



dW, 



jb 



1 



dr t 2W. 



J fc fc^j ° Tl I 2W jb Or, <Pb\ r 3k), 



a 



<f>b(rjk), 



if j = i 
if j ^ i 



(B4) 



The pairwise terms are given in Eq||, the derivatives are straightforward, noting that an atom does not interact with 
itself so V(ru) = 
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On combining the above equations, the force may then be written: 



E 



n ji ( n 3i A dMrj*) n j2 ( n j2 A dfojrji) 



4Wj i V N- 



4^-2^2 



9r; 



9ri <9r; 
",j /",. ,\ \ - d<j>i(r ik ) n j2 ( n j2 



d(f> 2 (r ik ) 
dri 



(B5) 



+ %i E 



aVi(r <fc ) 
9ri 



+ ^2 



dV 2 (r ik ) 
dri 



To simplify the appearance of this equation we define 



W ib 



n ib (N b - n ib ) 
±N b W ih 



(B6) 



and the force becomes 



+E 



dri 



dri 



-n j2 - 



dV 2 (r M ) 



dri 



, ^ 0&(r ifc ) avi(r ifc ) ay 2 (ri fe ) 

+ ^*2 — ?c riji — ^- nj 2 



dri 



dri 



(B7) 



The summations in the second set of square brackets may be reindexed in terms of j's and included in the first set of 
square brackets to give 



9ri 



+ W j2 - 



dr-; 



-fiji 



<9r. 



-n j2 - 



+ wii + _ njl d _Ypl _ n . 2 



<9r 4 

^ (r«) 



(B8) 



<9r 4 9ri dri dri 

Since both </> and V are pair potentials <fi b (rij) = 4> b (rji) and H^ij) = V b (rji). The derivatives of (f> are 

where the <f/ b are scalar quantities and 4>' b ( r ij) — 4>' b ( r ji)- Similarly, the derivatives of V are 



where again the V b ' are scalar quantities. 
Finally, the force on atom i is then 



S = E (W*+W jl )ci>' 1 (r ij ) + 














(B9) 


- +n jl )V{(r ji ) 


- {n i2 + n j2 )V 2 (rji) 







